home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / CollocWt.cc < prev    next >
C/C++ Source or Header  |  1996-03-03  |  4KB  |  195 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <iostream.h>
  32.  
  33. #include "CollocWt.h"
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36.  
  37. extern "C"
  38. {
  39.   int F77_FCN (jcobi, JCOBI) (int&, int&, int&, int&, double&,
  40.                   double&, double*, double*, double*,
  41.                   double*);
  42.  
  43.   int F77_FCN (dfopr, DFOPR) (int&, int&, int&, int&, int&, int&,
  44.                   double*, double*, double*, double*,
  45.                   double*);
  46. }
  47.  
  48. // Error handling.
  49.  
  50. void
  51. CollocWt::error (const char* msg)
  52. {
  53.   (*current_liboctave_error_handler) ("fatal CollocWt error: %s", msg);
  54. }
  55.  
  56. CollocWt&
  57. CollocWt::set_left (double val)
  58. {
  59.   if (val >= rb)
  60.     {
  61.       error ("left bound greater than right bound");
  62.       return *this;
  63.     }
  64.  
  65.   lb = val;
  66.   initialized = 0;
  67.   return *this;
  68. }
  69.  
  70. CollocWt&
  71. CollocWt::set_right (double val)
  72. {
  73.   if (val <= lb)
  74.     {
  75.       error ("right bound less than left bound");
  76.       return *this;
  77.     }
  78.  
  79.   rb = val;
  80.   initialized = 0;
  81.   return *this;
  82. }
  83.  
  84. void
  85. CollocWt::init (void)
  86. {
  87.   // Check for possible errors.
  88.  
  89.   double wid = rb - lb;
  90.   if (wid <= 0.0)
  91.     {
  92.       error ("width less than or equal to zero");
  93.       return;
  94.     }
  95.  
  96.   int nt = n + inc_left + inc_right;
  97.  
  98.   if (nt < 0)
  99.     {
  100.       error ("total number of collocation points less than zero");
  101.       return;
  102.     }
  103.   else if (nt == 0)
  104.     return;
  105.  
  106.   Array<double> dif1 (nt);
  107.   double *pdif1 = dif1.fortran_vec ();
  108.  
  109.   Array<double> dif2 (nt);
  110.   double *pdif2 = dif2.fortran_vec ();
  111.  
  112.   Array<double> dif3 (nt);
  113.   double *pdif3 = dif3.fortran_vec ();
  114.  
  115.   Array<double> vect (nt);
  116.   double *pvect = vect.fortran_vec ();
  117.  
  118.   r.resize (nt);
  119.   q.resize (nt);
  120.   A.resize (nt, nt);
  121.   B.resize (nt, nt);
  122.  
  123.   double *pr = r.fortran_vec ();
  124.  
  125.   // Compute roots.
  126.  
  127.   F77_FCN (jcobi, JCOBI) (nt, n, inc_left, inc_right, Alpha, Beta,
  128.               pdif1, pdif2, pdif3, pr);
  129.  
  130.   int id;
  131.  
  132.   // First derivative weights.
  133.  
  134.   id = 1;
  135.   for (int i = 1; i <= nt; i++)
  136.     {
  137.       F77_FCN (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1,
  138.                   pdif2, pdif3, pr, pvect); 
  139.  
  140.       for (int j = 0; j < nt; j++)
  141.     A (i-1, j) = vect.elem (j);
  142.     }
  143.  
  144.   // Second derivative weights.
  145.  
  146.   id = 2;
  147.   for (int i = 1; i <= nt; i++)
  148.     {
  149.       F77_FCN (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1,
  150.                   pdif2, pdif3, pr, pvect); 
  151.  
  152.       for (int j = 0; j < nt; j++)
  153.     B (i-1, j) = vect.elem (j);
  154.     }
  155.  
  156.   // Gaussian quadrature weights.
  157.  
  158.   id = 3;
  159.   double *pq = q.fortran_vec ();
  160.   F77_FCN (dfopr, DFOPR) (nt, n, inc_left, inc_right, id, id, pdif1,
  161.               pdif2, pdif3, pr, pq);
  162.  
  163.   initialized = 1;
  164. }
  165.  
  166. ostream&
  167. operator << (ostream& os, const CollocWt& a)
  168. {
  169.   if (a.left_included ())
  170.     os << "left  boundary is included\n";
  171.   else
  172.     os << "left  boundary is not included\n";
  173.  
  174.   if (a.right_included ())
  175.     os << "right boundary is included\n";
  176.   else
  177.     os << "right boundary is not included\n";
  178.  
  179.   os << "\n";
  180.  
  181.   os << a.Alpha << " " << a.Beta << "\n\n"
  182.      << a.r << "\n\n"
  183.      << a.q << "\n\n"
  184.      << a.A << "\n"
  185.      << a.B << "\n";
  186.  
  187.   return os;
  188. }
  189.  
  190. /*
  191. ;;; Local Variables: ***
  192. ;;; mode: C++ ***
  193. ;;; End: ***
  194. */
  195.